#include "NormalQApprox.h"

#include <math.h>
#include <errno.h>


NormalQApprox::NormalQApprox(double mu, double sigma) : LOW(0.02425), HIGH( 0.97575), _mu(mu), _sigma(sigma)
{

	/* Coefficients in rational approximations. */
	static const double arrA[] = {
		-3.969683028665376e+01,
		2.209460984245205e+02,
		-2.759285104469687e+02,
		1.383577518672690e+02,
		-3.066479806614716e+01,
		2.506628277459239e+00
	};

	static const double arrB[] =
	{
		-5.447609879822406e+01,
		1.615858368580409e+02,
		-1.556989798598866e+02,
		6.680131188771972e+01,
		-1.328068155288572e+01
	};

	static const double arrC[] =
	{
		-7.784894002430293e-03,
		-3.223964580411365e-01,
		-2.400758277161838e+00,
		-2.549732539343734e+00,
		4.374664141464968e+00,
		2.938163982698783e+00
	};

	static const double arrD[] =
	{
		7.784695709041462e-03,
		3.224671290700398e-01,
		2.445134137142996e+00,
		3.754408661907416e+00
	};

	a.assign(arrA, arrA + sizeof(arrA) / sizeof(arrA[0]) );
	b.assign(arrB, arrB + sizeof(arrB) / sizeof(arrB[0]) );
	c.assign(arrC, arrC + sizeof(arrC) / sizeof(arrC[0]) );
	d.assign(arrD, arrD + sizeof(arrD) / sizeof(arrD[0]) );

}

NormalQApprox::~NormalQApprox(void)
{
}

/*
	* The inverse standard normal distribution.
	*
	*   Author:      Peter John Acklam <pjacklam@online.no>
	*   URL:         http://home.online.no/~pjacklam
*/

double NormalQApprox::operator()( double p )
{
	double xApprox;
	double q, r;

	errno = 0;

	if (p < 0 || p > 1)
	{
		errno = EDOM;
		return 0.0;
	}
	else if (p == 0)
	{
		errno = ERANGE;
		return -HUGE_VAL /* minus "infinity" */;
	}
	else if (p == 1)
	{
		errno = ERANGE;
		return HUGE_VAL /* "infinity" */;
	}
	else if (p < LOW)
	{
		/* Rational approximation for lower region */
		q = sqrt(-2*log(p));
		xApprox = (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
			((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1);
	}
	else if (p > HIGH)
	{
		/* Rational approximation for upper region */
		q  = sqrt(-2*log(1-p));
		xApprox = -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
			((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1);
	}
	else
	{
		/* Rational approximation for central region */
		q = p - 0.5;
		r = q*q;
		xApprox = (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q /
			(((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1);
	}

	return _mu + _sigma*xApprox;
}
 
